function ErrorCalculatorFebST
format long e

gauge1=load('febx50.out');
gauge2=load('febx100.out');
gauge3=load('febx150.out');
gauge4=load('febx200.out');
gauge5=load('febx250.out');
gauge6=load('febx300.out');
gauge7=load('febx350.out');
gauge8=load('febx400.out');
gauge9=load('febx450.out');
gauge10=load('febx500.out');
gauge11=load('febx550.out');
gauge12=load('febx600.out');
gauge13=load('febx650.out');

datax50=[gauge1(:,1) gauge1(:,2)];
datax100=[gauge2(:,1) gauge2(:,2)];
datax150=[gauge3(:,1) gauge3(:,2)];
datax200=[gauge4(:,1) gauge4(:,2)];
datax250=[gauge5(:,1) gauge5(:,2)];
datax300=[gauge6(:,1) gauge6(:,2)];
datax350=[gauge7(:,1) gauge7(:,2)];
datax400=[gauge8(:,1) gauge8(:,2)];
datax450=[gauge9(:,1) gauge9(:,2)];
datax500=[gauge10(:,1) gauge10(:,2)];
datax550=[gauge11(:,1) gauge11(:,2)];
datax600=[gauge12(:,1) gauge12(:,2)];
datax650=[gauge13(:,1) gauge13(:,2)];

NNN=length(gauge1);
fftx50=fft(datax50(:,2))/NNN;
fftx100=fft(datax100(:,2))/NNN;
fftx150=fft(datax150(:,2))/NNN;
fftx200=fft(datax200(:,2))/NNN;
fftx250=fft(datax250(:,2))/NNN;
fftx300=fft(datax300(:,2))/NNN;
fftx350=fft(datax350(:,2))/NNN;
fftx400=fft(datax400(:,2))/NNN;
fftx450=fft(datax450(:,2))/NNN;
fftx500=fft(datax500(:,2))/NNN;
fftx550=fft(datax550(:,2))/NNN;
fftx600=fft(datax600(:,2))/NNN;
fftx650=fft(datax650(:,2))/NNN;
fftxAll=[fftx50,fftx100,fftx150,fftx200,fftx250,fftx300,fftx350,fftx400,fftx450,fftx500,fftx550,fftx600,fftx650]';

ZeroModeLoc=101;
NumModes=50;
% The following k0 value includes surface tension.
k0=0.440938548760477444;
epsilon=2*k0*abs(fftx50(ZeroModeLoc));

ExptGauges=abs(fftxAll(:,ZeroModeLoc-NumModes:ZeroModeLoc+NumModes));
[MMM,JJJ]=size(ExptGauges);

NLSSTAlldata=epsilon/k0*load('NLSSTAlldata.out');
dNLSSTAlldata=epsilon/k0*load('dNLSSTAlldata.out');
DystheSTAlldata=epsilon/k0*load('DystheSTAlldata.out');
vDystheSTAlldata=epsilon/k0*load('vDystheSTAlldata.out');
dGTSTAlldata=epsilon/k0*load('dGTSTAlldata.out');
NLSSTGauges=NLSSTAlldata(1:500:6001,:);
dNLSSTGauges=dNLSSTAlldata(1:500:6001,:);
DystheSTGauges=DystheSTAlldata(1:500:6001,:);
vDystheSTGauges=vDystheSTAlldata(1:500:6001,:);
dGTSTGauges=dGTSTAlldata(1:500:6001,:);

HannahNLSSTError=100/JJJ/(MMM-1)*sum(sum((ExptGauges'-NLSSTGauges').^2)./sum((abs(ExptGauges')).^2))
HannahdNLSSTError=100/JJJ/(MMM-1)*sum(sum((ExptGauges'-dNLSSTGauges').^2)./sum((abs(ExptGauges')).^2))
HannahDystheSTError=100/JJJ/(MMM-1)*sum(sum((ExptGauges'-DystheSTGauges').^2)./sum((abs(ExptGauges')).^2))
HannahvDystheSTError=100/JJJ/(MMM-1)*sum(sum((ExptGauges'-vDystheSTGauges').^2)./sum((abs(ExptGauges')).^2))
HannahdGTSTError=100/JJJ/(MMM-1)*sum(sum((ExptGauges'-dGTSTGauges').^2)./sum((abs(ExptGauges')).^2))

u1p=3;
u2p=7;
u3p=10;
l1m=3;
l2m=7;
l3m=10;

a0s=abs(fftxAll(:,ZeroModeLoc));
an1s=abs(fftxAll(:,ZeroModeLoc-l1m));
a1s=abs(fftxAll(:,ZeroModeLoc+u1p));
an2s=abs(fftxAll(:,ZeroModeLoc-l2m));
a2s=abs(fftxAll(:,ZeroModeLoc+u2p));
an3s=abs(fftxAll(:,ZeroModeLoc-l3m));
a3s=abs(fftxAll(:,ZeroModeLoc+u3p));

carrier=(JJJ+1)/2;
u1=carrier+u1p;
u2=carrier+u2p;
u3=carrier+u3p;
l1=carrier-l1m;
l2=carrier-l2m;
l3=carrier-l3m;

NLSSTcarrier=NLSSTAlldata(1:100:end,carrier);
dNLSSTcarrier=dNLSSTAlldata(1:100:end,carrier);
DystheSTcarrier=DystheSTAlldata(1:100:end,carrier);
vDystheSTcarrier=vDystheSTAlldata(1:100:end,carrier);
dGTSTcarrier=dGTSTAlldata(1:100:end,carrier);
NLSu1=NLSSTAlldata(1:100:end,u1);
dNLSu1=dNLSSTAlldata(1:100:end,u1);
Dystheu1=DystheSTAlldata(1:100:end,u1);
vDystheu1=vDystheSTAlldata(1:100:end,u1);
dGTSTu1=dGTSTAlldata(1:100:end,u1);
NLSu2=NLSSTAlldata(1:100:end,u2);
dNLSu2=dNLSSTAlldata(1:100:end,u2);
Dystheu2=DystheSTAlldata(1:100:end,u2);
vDystheu2=vDystheSTAlldata(1:100:end,u2);
dGTSTu2=dGTSTAlldata(1:100:end,u2);
NLSu3=NLSSTAlldata(1:100:end,u3);
dNLSu3=dNLSSTAlldata(1:100:end,u3);
Dystheu3=DystheSTAlldata(1:100:end,u3);
vDystheu3=vDystheSTAlldata(1:100:end,u3);
dGTSTu3=dGTSTAlldata(1:100:end,u3);

NLSl1=NLSSTAlldata(1:100:end,l1);
dNLSl1=dNLSSTAlldata(1:100:end,l1);
Dysthel1=DystheSTAlldata(1:100:end,l1);
vDysthel1=vDystheSTAlldata(1:100:end,l1);
dGTSTl1=dGTSTAlldata(1:100:end,l1);
NLSl2=NLSSTAlldata(1:100:end,l2);
dNLSl2=dNLSSTAlldata(1:100:end,l2);
Dysthel2=DystheSTAlldata(1:100:end,l2);
vDysthel2=vDystheSTAlldata(1:100:end,l2);
dGTSTl2=dGTSTAlldata(1:100:end,l2);
NLSl3=NLSSTAlldata(1:100:end,l3);
dNLSl3=dNLSSTAlldata(1:100:end,l3);
Dysthel3=DystheSTAlldata(1:100:end,l3);
vDysthel3=vDystheSTAlldata(1:100:end,l3);
dGTSTl3=dGTSTAlldata(1:100:end,l3);

positions=[0,50,100,150,200,250,300,350,400,450,500,550,600]';
x=0:1/10:6000/10;
x=x(1:100:end);

figure(6)
subplot(4,2,[1,2])
hold off
scatter(positions,a0s,'filled','k')
hold on
plot(x,NLSSTcarrier,'k')
plot(x,DystheSTcarrier,'k--')
l=legend('data','NLS','Dysthe','Location','SouthWest','Orientation','horizontal');
l.NumColumns=2;
ylabel('$|a_0|$','Interpreter','Latex')
title('Expt C First Harmonic')
axis tight
subplot(4,2,3)
hold off
scatter(positions,an1s,'filled','k')
hold on
plot(x,NLSl1,'k')
plot(x,Dysthel1,'k--')
ylabel('$|a_{-3}|$','Interpreter','Latex')
xlim([0 x(end)])
ylim([0 0.05])
subplot(4,2,4)
hold off
scatter(positions,a1s,'filled','k')
hold on
plot(x,NLSu1,'k')
plot(x,Dystheu1,'k--')
ylabel('$|a_3|$','Interpreter','Latex')
xlim([0 x(end)])
ylim([0 0.05])
subplot(4,2,5)
hold off
scatter(positions,an2s,'filled','k')
hold on
plot(x,NLSl2,'k')
plot(x,Dysthel2,'k--')
ylabel('$|a_{-7}|$','Interpreter','Latex')
xlim([0 x(end)])
ylim([0 0.03])
subplot(4,2,6)
hold off
scatter(positions,a2s,'filled','k')
hold on
plot(x,NLSu2,'k')
plot(x,Dystheu2,'k--')
ylabel('$|a_7|$','Interpreter','Latex')
xlim([0 x(end)])
ylim([0 0.03])
subplot(4,2,7)
hold off
scatter(positions,an3s,'filled','k')
hold on
plot(x,NLSl3,'k')
plot(x,Dysthel3,'k--')
xlabel('$x$ (cm)','Interpreter','Latex')
ylabel('$|a_{-10}|$','Interpreter','Latex')
xlim([0 x(end)])
ylim([0 0.02])
subplot(4,2,8)
hold off
scatter(positions,a3s,'filled','k')
hold on
plot(x,NLSu3,'k')
plot(x,Dystheu3,'k--')
xlabel('$x$ (cm)','Interpreter','Latex')
ylabel('$|a_{10}|$','Interpreter','Latex')
xlim([0 x(end)])
ylim([0 0.02])



figure(66)
subplot(4,2,[1,2])
hold off
scatter(positions,a0s,'filled','k')
hold on
plot(x,dNLSSTcarrier,'k')
plot(x,vDystheSTcarrier,'k--')
plot(x,dGTSTcarrier,'k-.')
l=legend('data','dNLS','vDysthe','dGT','Location','SouthWest','Orientation','horizontal');
l.NumColumns=2;
ylabel('$|a_0|$','Interpreter','Latex')
title('Expt C First Harmonic')
axis tight
subplot(4,2,3)
hold off
scatter(positions,an1s,'filled','k')
hold on
plot(x,dNLSl1,'k')
plot(x,vDysthel1,'k--')
plot(x,dGTSTl1,'k-.')
ylabel('$|a_{-3}|$','Interpreter','Latex')
xlim([0 x(end)])
ylim([0 0.05])
subplot(4,2,4)
hold off
scatter(positions,a1s,'filled','k')
hold on
plot(x,dNLSu1,'k')
plot(x,vDystheu1,'k--')
plot(x,dGTSTu1,'k-.')
ylabel('$|a_3|$','Interpreter','Latex')
xlim([0 x(end)])
ylim([0 0.05])
subplot(4,2,5)
hold off
scatter(positions,an2s,'filled','k')
hold on
plot(x,dNLSl2,'k')
plot(x,vDysthel2,'k--')
plot(x,dGTSTl2,'k-.')
ylabel('$|a_{-7}|$','Interpreter','Latex')
xlim([0 x(end)])
ylim([0 0.03])
subplot(4,2,6)
hold off
scatter(positions,a2s,'filled','k')
hold on
plot(x,dNLSu2,'k')
plot(x,vDystheu2,'k--')
plot(x,dGTSTu2,'k-.')
ylabel('$|a_7|$','Interpreter','Latex')
xlim([0 x(end)])
ylim([0 0.03])
subplot(4,2,7)
hold off
scatter(positions,an3s,'filled','k')
hold on
plot(x,dNLSl3,'k')
plot(x,vDysthel3,'k--')
plot(x,dGTSTl3,'k-.')
xlabel('$x$ (cm)','Interpreter','Latex')
ylabel('$|a_{-10}|$','Interpreter','Latex')
xlim([0 x(end)])
ylim([0 0.02])
subplot(4,2,8)
hold off
scatter(positions,a3s,'filled','k')
hold on
plot(x,dNLSu3,'k')
plot(x,vDystheu3,'k--')
plot(x,dGTSTu3,'k-.')
xlabel('$x$ (cm)','Interpreter','Latex')
ylabel('$|a_{10}|$','Interpreter','Latex')
xlim([0 x(end)])
ylim([0 0.02])
